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1.0  Introduction 


VLSI  (Very  Largo  Scale  Integration)  technology  has  been  developed  to  the 
point  where  high  speed  floating  point  processors  may  bo  concatenated  to  form 
compact  supercomputers  with  far  greater  throughput  than  uniprocessor  machines. 
Thus,  there  currently  is  considerable  interest  within  the  signal  processing 
community  in  the  development  of  parallel  versions  of  many  conventional 
algorithms  such  as  convolution,  matrix  multiplication  and  factorization.  MTI 
has  collaborated  with  Dr.  G.J.  Bierman  to  develop  a  parallel  form  of  the 
Kalman  filter  that  has  several  very  unique  and  important  features.  We  believe 
that  utilization  of  these  features  will  result  in  the  design  of  an  integrated 
test  range  tracking  system  that  exhibits  much  improved  performance  over  the 
existing  "independent"  approach  to  tracking  at  the  White  Sands  Missile  Range 
( WSMR ) . 

Specifically,  our  Decentralized  Square  Root  Information  Filter  (DSRIF) 

[1]  allows  each  group  of  measurement  variables,  the  process  noise  statistics 
and  the  prior  information  about  the  initial  state  to  be  processed  in  separate 
but  locally  optimal  filters.  Globally  optimal  state  estimates  and  estimate 
error  covariances  may  then  be  computed  by  combining  local  filter  outputs  on 
demand.  This  will  allow  thr  analyst  to  identify  the  contribution  of  each 
measurement  group,  the  process  noise  and  the  prior  information  about  the 
initial  state  ta  the  global  state  estimate  and  estimate  error  covariance 
without  additional  computation. 

Furthermore,  the  process  noise  and  prior  information  may  bo  distributed 
amongst  the  data  processing  filters  in  order  to  improve  upon  the  fault 
tolerant  characteristics  of  the  nominal  olgorithm  when  real-time  signal 
processing  is  an  issue.  In  this  case,  the  estimates  and  covariances  should 
gracefully  degrade  from  global  optimality  as  local  processors  foil.  Thirdly, 
the  algorithm  is  based  upon  numerically  reliable  matrix  factorization  methods 
which,  unlike  the  Conventional  Kalman  Filter  (CKF),  will  never  fail. 

The  objective  of  Phase  I  research  was  to  validate  the  DSRIF  equations  by 
testing  its  ability  to  track  a  predetermined  ballistic  trajectory  when 
perturbed  by  white  Gaussian  noise.  The  stote  estimates  and  error  covariances 
obtained  were  found  to  be  identical  (when  printed  to  10  significant  digits) 
with  those  of  a  SRIF  implemented  in  centralized  form  with  all  calculations 
performed  in  double  precision  arithmetic.  Furthermore,  an  extended  version  of 
the  DSRIF  (E-DSRIF)  was  derived  and  successfully  used  to  track  real  Multiple 
Rocket  Launch  System  data  obtained  from  the  WSMR. 

The  remaining  parts  of  this  Phase  I  report  ore  orgonized  as  follows.  In 
the  next  section  we  formulate  the  decentralized  estimation  problem  and  define 
a  necessary  and  sufficient  condition  for  recovering  globally  optimal  state 
estimates  from  locally  optimol  ones.  Then,  in  section  2.1  wo  derive  the  DSRIF 
and  show  how  the  effects  of  prior  and  process  noise  may  bo  distributed.  In 
section  2.2,  the  ability  of  the  DSRIF  to  accurately  track  the  position  and 
velocity  of  a  ballistic  object  is  tested  via  a  computer  simulation.  Finally, 
in  section  2.3,  the  E-DSRIF  is  used  to  track  a  maneuvering  vehicle  but  using  a 
polynomial  process  model  for  the  target  dynamics. 
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2 . 0  Results  of  the  Phase  I  Work 


The  trajectory  estimation  problem  is  a  problem  of  nonlinear  estimation. 
Assume  that  the  state  of  a  target  evolves  in  time  according  to  the  equation 

x(t)  -  f(x{t),u(t))  +  w(t)  (1) 

where  u(t)  is  the  target's  nominal  control  vector  and  w(t)  is  a  zero  mean 
white  noise  process  with  spectral  density  Q(t).  Values  for  the  latter  are 
selected  in  order  to  compensate  for  errors  in  the  model  which  may  originate 
from  unknown  perturbations  to  the  nominal  control  (such  as  shear  winds,  wind 
gusting  etc.)  as  well  as  from  uncertainties  in  the  aerodynamics  of  the  target 
vehicle.  The  corresponding  discrete  measurement  vector  for  the  i^^  sensor  is 
given  by 


Vk  =  ^^(x,^)  +  vj  (2) 

where  v^  is  a  zero  mean  white  noise  sequence  with  covariance  Rj.  The  problem 
is  to  estimate  the  target  states  x^  based  upon  all  of  the  past  measurements  y^^ 
where  {l<i<M  and  1<l<k}.  The  state  vector  X(^  contains  the  target  position  and 
velocity,  biases  which  account  for  the  displacement  and  orientation  of  the 
sensor  or  "local"  coordinate  systems  (LCSs)  w.r.t.  the  global  coordinate 
system  (GCS),  and  acceleration  when  the  target  vehicle  is  maneuvering. 

Unfortunately,  the  optimal  nonlinear  estimator  (conditional  moan)  cannot 
be  realized  with  a  finite-dimensional  implementation  and  consequently,  all 
practical  nonlinear  filters  must  be  suboptimal.  The  usual  suboptimal  solution 
is  the  CKF  when  the  nominal  trajectory  is  known  o  priori,  the  Extended  Kalman 
Filter  when  the  nominal  trajectory  is  unavailable,  and  Higher  Order  Filters 
(such  as  the  second  order  filter  [4],  the  single-stage  iterative  filter  [5], 
and  the  Gaussian  sum  filter  [6],  among  others  [7])  when  even  greater  accuracy 
is  desired.  The  tradeoff  here  is  performance  versus  real-time  computational 
requirement. 

Thus,  wo  see  that  the  ability  of  a  sensor  group  to  accurately  record  the 
motion  of  one  or  more  airborne  targets  is  a  function  of  the  individual  target 
and  sensor  dynamical  models  as  well  as  the  particular  algorithm  used  to 
combine  raw  data  and  produce  track  estimates.  Another  issue  is  intersensor 
communication.  Sensors  which  operate  independently  from  one  another  will 
exhibit  larger  estimate  errors  than  ones  which  communicate  with  other  members 
of  the  network. 

Let  the  global  discrete  time  linear  system 

Xk  •  ^k-1  ’^k-l  ''^k-1 

where 

w^  -  N(0,0^)  (4) 

Xq  -  N(0,Po(-))  (5) 

be  the  model  for  the  target,  and  the  global  measurement  model 


+ 


(6) 


— 

Yk 

a 

■ 

r 

< 

L 

o 

J 

where  VjJ  through  vJJ  ore  uncorrelated  random  vectors  and 

vj  .  N(O.Rj)  (7) 

be  the  model  for  the  sensors."  The  problem  is  to  calculate  the  globally 
optimal  (minimum  mean  square  error)  estimate  of  Xj^  and  its  associated  estimate 
error  covariance  matrix  when  through  yJJ  are  processed  separately  by 
locally  optimal  estimators  1  through  M  correspondingly.  The  local  dynamical 
models  are 

’‘k  ■  ^k-1  ’‘k-1  '^k-1 

where 

w^  -  N{O.Qj) 
xj  -  N{0.pj(-)) 

and  the  local  measurement  models  are 

y^  -  Hj  xj  +  v^ 

where  vj  satisfies  (7).  Notice  that  the  local  states  may  be  physically 
different  from  the  global  states.  Wilsky  et.al.  [3]  have  recently  shown  that 
a  necessary  and  sufficient  condition  for  our  being  able  to  recover  globally 
optimal  state  estimates  from  locally  optimal  ones  is  that 

cj  .  4  4  (12) 

and  for  the  tracking  application,  we  expect  that 

x^  -  Xk  (13) 

where  is  a  matrix  which  results  in  the  correct  partitioning  of  global 
states  to  the  subsystem  filters.  Alternatively,  equation  (12)  allows  us  to 
define  the  local  state  vector  os  d^  plus  the  local  state  vector  defined  in 
section  2.2  (both  in  local  coordinates)  provided  that  wj  =  (T^)^'". 


(8) 

(9) 

(10) 

(11) 


If 

Lower  case  variables  are  vector  quantities  while  upper  case  generally 
corresponds  to  matrices  of  appropriate  dimension.  Also  w^,  xq,  and  V|J  through 
vjJ  are  uncorrelated  with  each  other  for  all  k  and  N(a,B)  signifies  an  a  mean 
white  Gaussian  process  with  covariance  matrix  B. 
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The  derivation  of  this  algorithm  may  be  found  in  our  recent  publication 
[1]  with  the  exception  of  some  new  ideas  on  the  distribution  of  prior 
information  about  the  process  (plont)  noise  and  initial  state  amongst  the 
local  processors. 


Let  =  I.  The  goal  is  to  find  the  sequence  xq . xjg  that  minimizes 

the  least  squares  likelihood  performance  functional  [2,  pg.42] 


. °  II  **0  “  ^0^“^  I  l‘ 


M  N 

Z  Z 
i=1  k=0 


(Ri)-^Hi  X;,  - 


+  E  I  I  z^(  k  )  -  R^^{  k )  w^ 
k-0 


where  the  a  priori  estimate  Rg  (-)  Zq(-)  has  covariance 


Pfl(-) 


(-)  Rfl  (-)  , 


-1  -T 

Qk  -  Rw  (k)  Rw  (k)  . 


Zv^(k)  «  R^(k)  times  the  a  priori  expectation  of  w^,  ond  the  pair  zj,  hj 
correspond  to  normalized  measurement  equations  i.e.,  Pj  »  I.  Decentralized 
processing  is  achieved  by  distributing  the  minimization  of  the  performance 
criterion  amongst  the  local  filters  ond  global  merging  equations  that  follow. 
The  best  distribution  for  target  tracking  is  probably  to  minimize 
[ I  Hj  xk  -  zj  ! 1^  in  each  of  the  M  local  filters  and  minimize  the  remaining 
two  terms  in  (14)  within  the  central  (merge)  processor.  However,  this  point 
should  bo  explored  in  detail  in  Phase  II  research.  Thus,  both  data  types  may 
be  processed  using  the  standard  SRIP  mechanization: 


Measurement  Update 


Ri(-) 

zi(-) 

Rk(+)  zi(+) 

(Ri)-=^  Hi 

(Ri)-%  yi 

r 

o 

o 

1 _ 

Time  Update 


Rw(k) 


Zw(k) 


-Ri(+)  Ri(+)  V 


zi(  +  ) 


(i) 


R^Ck) 


Rwx(k) 


(i) 

(i) 


z^Ck) 


(i) 

(i) 


Rk+i(-) 


^k+l("5 


(18) 


where 


[  Rj(-)  zj(-)  ]  -  [  0  0  ] 


[  Rv^(k)  z^^(k)  ]  -  [  0  0  ] 


(19) 

(20) 


and  are  orthogonal  transformations  which  put  the  matrices  on  the  left 


hand  sides  of  (17)  and  (18)  into  upper  triangular  form.  They  may  be 
implicitly  computed  using  Householder  transformations  [2,  p. 60-64]. 


In  terms  of  the  local  filter  results  (17)  and  (18),  the  performance  can 
bo  rewritten  as 

N-1 

Jn(’'o . *  II  ''0  ■  11^  +  ^  II  -  Rw^*')  '"k  11^ 

k-0 


(1) 

1 

(1) 

r  — 1 

(1) 

Rw(K) 

Wk  + 

fCx^K) 

^k+1  - 

z*(k) 

(M) 

(M) 

(M) 

r 

5r 

L 

fCx(K) 

z:(k) 

M  M 
+  E  E 
i=1  k-0 


l|2  t  II 


r~  n 

(1) 

(1) 

Rn(+) 

Zn(  +  ) 

• 

’‘N  - 

(M) 

(M) 

Rn(  +  ) 

Zn(  +  ) 

(21) 


Applying  an  orthogonal  transformation  to  the  first  three  terms  in  (21), 
results  in  the  following  recursive  equotion  for  combining  local  smoothing 
coefficients  with  process  noise  and  prior  on  Xg: 
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Z  (k-1)  J 


I  -H*(k-1)  ij, 


fC(k)  RCx^''> 

0  H"(k)  z^’Ck) 

0  0# 


(22 


[  H*(-1)  z^-l)  ]  =  C  Ro(-)  zo(-)  ]  (23 

and  H*(k)  is  upper  triangular.  To  obtain  the  globally  optimal  information 
vector  Z(,(+)  and  square  root  information  matrix  R^(+),  wo  solve  the  following 
equations  using  H*(k)  and  z*(k)  from  (22); 


(1)  (1) 
Rk(  +  ) 


(M)  (M) 

R^(  +  ) 

H*(k-1)  z*(k-1) 


R^(+)  z^(+) 

0  # 


(24 


where  is  an  orthogonal  transformation  which  puts  the  loft  hand  side  of 
(24)  in  upper  triangular  form.  Globally  optimal  filter  estimates  and 
covariances  are  then  given  by 


-1 

Xk(+)  -  R^  (+)  z^(+)  (25 

-1  -tr 

P^(  +  )  -  R,^  (  +  )  R^  (  +  )  (26 

When  the  a  priori  information  about  the  initial  condition  and  process  noise 
models  are  adjusted,  it  is  only  necessary  to  rerun  (22)  and  (24)  without 
reprocessing  any  measurements. 


To  summarize,  each  data-typo  is  processed  by  a  local  SRIF 


^(i) 

(17)-(20)  which  generates  a  set  of  smoothing  coefficients  R^  , 

(i)  (i)  (i> 

R^(+)  ,  z*(+)  as  well  as  a  square  root  information  motrix  R(+)  and 

(i) 

information  vector  z(+)  The  central  or  merge  processor  consists  of  three 

separote  processors  which  operate  in  porollel.  The  first  mechonizes  (22), (23) 
which  combines  the  local  smoothing  coefficients  with  the  effects  of  process 
noise  and  prior  information  about  the  initial  state.  The  second  mechanizes 
(24)  which  merges  the  local  square  root  information  matrices  and  vectors  with 
output  from  the  first,  but  only  upon  demand  by  the  third.  The  third  produces 
estimates  and  covariances  whenever  desired  by  back-solving  (25)  and  (26) 
respectively,  noting  that  (25), (26)  require  output  from  the  second  processor. 
An  important  observation  is  that  there  is  no  feedbock  of  information  from  the 
merge  processor  to  the  local  filters. 

Actually,  a  family  of  DSRIFs  exist.  Each  member  corresponds  to  a 
different  distribution  of  process  noise  and  prior  information  about  the 
initial  state  amongst  the  local  processors  and  the  first  merge  processor 
(22), (23).  To  see  this,  realize  that  the  prior  and  process  noise  terms  in 
(14)  may  be  written  as 


M 

(i) 

(i) 

Zw(k)  -  R^(k)  w,, 

11^  ■  E  1 

i»1 

1  ZwO') 

-  Ry,(  k )  w,^ 

M 

(i) 

(i) 

11^  ■  I  1 

!  Ro(-) 

X 

o 

1 

N 

O 

1 

i»1 


(i)  (i) 

.zo(-) 


where  the  dimension  of  w  and  x  are  both  n.  The  vectors  z^(k) 


(i)  (i) 

matrices  R^(k)  .RgC”)  af®  chosen  so  that  the  local  norms  on  the  right 
hand  side  of  (27), (28)  remain  invariant  with  respect  to  orthogonal 
transformations.  Each  norm  may  be  incorporated  into  a  local  SRIF  or  the  first 
merge  processor  and  minimized  there. 

2.2  Validation  of  the  DSRIF  for  a  Simulated  Ballistic 
Traiectorv  Over  the  White  Sands  Missile  Ranae 


In  order  to  validate  our  decentralized  approach  to  solving  the  linear 
least  squares  estimation  problem,  a  high  fidelity  simulation  of  a  multisensor 
network,  tracking  a  ballistic  trajectory  over  the  WSMR,  was  encoded  in  Fortran 
'77  and  executed  on  an  IBM  (clone)  Model  AT  desktop  computer  (640K  ram,  20 
Mbyte  hard  disk,  Intel  80287  math  coprocessor).  Initial  conditions  for  the 
nominal  trajectory  were  calculated  using  a  "flat  earth"  or  constant 
gravitational  acceleration  model  which  neglects  the  earth’s  rotation  and 
assumes  that  the  target  vehicle  is  a  point  mass. 

The  initial  position  is  launch  complex  #32  (E4,N2)  and  the  desired 
terminal  position  is  the  GAM  83  target  (E13,N85).  This  corresponds  to  a 
flight  path  of 

((13-4)2  +  (85-2)2)’^  -  83.49  miles  (29) 

when  projected  onto  the  "Range  Area  -  General  Road  Mop  (RAGRM)"  which  we  use 
to  define  the  GCS.  That  is  x^,x2,x^  is  a  right  handed  coordinate  system  where 
the  vector  cross  product  of  x^  with  x2  -  x^  and  the  vectors  [x^  0  0]  and  [0  x2 
0]  point  east  along  latitude  32.380  degrees  and  north  along  longitude  106.481 
degrees  respectively.  Then  [0  0  x^]  is  collinear  with  the  rodial  vector  which 
points  outward  from  the  earth’s  center  and  passes  through  the  origin  of  the 
GCS.  A  launch  elevation  angle,  w.r.t.  the  x\x2  tangent  plane,  of  45  degrees 
was  chosen  in  order  to  maximize  the  projected  flight  path  for  a  given  amount 
of  energy.  The  corresponding  initial  positions  and  velocities  are  then 

xj  -  21 , 120.  feet 

Xq  >  10,560.  feet 

xg  »  0.  feet 

xg  -  ((32.  feet/sec2)(83 . 49  miles/2. )(5280  feet/mile))^  (30) 

-  2655.8  feet/sec 

Xq  "  (x§)(9./83.49)  -  286.29  feet/sec 
Xq  -  (x§)(83./83.49)  -  2640.21  feet/sec 
with  a  maximum  ensuing  altitude  of 

x5(tf^/2.)  -  (x§)2/(2.(32.  foet/sec2))  .  110,207.  feet  (20.9  miles) 

(31) 

where 
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165.99  seconds, 


(32) 


»  (83.49  miles)(5280  feet/tnile)/x®(0) 
the  flight  impact  time. 

Five  sensor  locations  which  cover  both  sides  of  the  projected  flight  path 
were  selected.  They  are  Zebra  (E3,N14),  Rhodes  Canyon  Range  Center  (E0,N53), 
A8RES  Radar  (W10,N99),  Oscuro  Range  Center  (E18,N77)  and  King  I  (E21,N34). 

Each  sensor  records  measurements  with  respect  to  its  own  LCS  so  that 
coordinate  transformations  to  the  GCS  were  derived  and  included  in  the 
observational  equations.  The  transformation  is 


x! 

r-  .  .-I 

x^’  ^ 

x2 

x3 

_  _ 

-  (Ti) 

xi,2 

xi.3 

(33) 


where 


(  ®Ti  )  (  ^Ti  )  (  ®T^  )  (  ®Ti  ) 


(34) 


and 


Sfi 


1  0  0 

0  cos(L+a^)  sin(L+a^) 

0  -sin(L+a^)  cos(L+a^) 


(35) 


Sfi 


cos  0 

0  .  1 

sin  p^  0 


-sin  P'* 
0 

cos  p^ 


(36) 


7Ti 


Sri 


0 

cos  L 
sin  L 


-sin  T 
0 


sin 

cos 

0 


0 

-sin  L 
cos  L 


(37) 


(38) 


noting  that  (T^)^*"  -  (T^)"^  since  T^  is  an  prthogonal  transformation.  A 
computational  savings  results  when  (®T^)(®T^ )(^T^ )  are  combined  using 
trigonometric  identities  for  the  cosine  and  sine  of  the  sum  of  two  angles 
along  with  the  small  angle  approximation  since  is  bounded  by  +  1  degrees 
over  the  length  of  the  WSMR. 


d^  Is  the  vector  from  the  GCS  to  the  ith  LCS  and  a^.p^.r^  are  the  three 
Euler  angles  describing  the  orientation  of  the  ith  LCS  w.r.t.  the  GCS. 
Specifically,  if  we  first  rotate  about  counterolockwise  by  an  (L+a^)  degree 
change  in  latitude  (aligning  y^  with  the  polar  axis),  and  then  rqtate 


counterclockwise  about  y^ 


y^  by  a  p^  degree  change  in  longitude,  (p^  is  bounded 


by  +  1/2  degree  over  the  width  of  the  WSMR)  then  rotate  clockwise  about  x3  by 
on  L  degree  change  in  latitude,  and  finally  rotate  about  by  degrees  to 
account  for  "tangent  plane  misalignment”,  the  LCS  will  coincide  with  the  GCS 
when  d^  »  0.  Is  equal  to  the  LCS  latitude  -  GCS  latitude,  Is  equal  to 
the  LCS  longitude  -  GCS  longitude. 


To  determine  a^,  given  L, 
following  two  equations: 


S  and  d^ 


(in  global  coordinates). 


solve  the 


ui.  1 1 
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m 
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x"' 

°  2 

rg  cos  L  cos  S 

■ 

-rg  cos  L  sin  S 

OX^ 

ECEF 

rg  sin  L 
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— 

~1 

— 

— 
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S 

-cos  S 
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1 
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S 

sin  S 

0 

0 

sin  L 
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L 

0 
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1 
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(40) 


(41) 


1.1 

1.2 

1.3 


ECEF 


gr^  cos  (L+af)  cos  (S+p^) 
gr^  cos  (L+a^)  sin  (S+p^) 
gr^  sin  (L+a^) 


(42) 


The  solution  is 


L+a^  »  tan~^ 


(gX^>‘')2  +  (gX^’2)2 


(L+a^)  s  32  degrees  (43) 


S+p^  »  cos' 


gr^  cos  (L+a^) 


(S+p^)  s  106  degrees 


(44) 


o 


X 


i.3 


where 


The  results  are  given  in  Table  1  below. 


I 
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Table  1:  Optical  Tracker  (ot)  and  Range  Radar  (rr)  Locations  o’-d 
Orientations  in  GCS  for  the  Ballistic  Target 

Detailed  equatians  that  describe  the  translational  motion  of  the  target 
were  developed.  The  equations  include  a  radial  gravitational  force  as  well  as 
centrifugal  and  Coriolis  forces  which  come  about  by  rotation  of  the  GCS  about 
the  polar  axis.  The  equations  are 


((x1)2  +  (x2)2  +  (x3+rg)2)3/2 
->ix2 

((x1)2  +  (x2)2  +  (x3+r^)2)3/2 

-p(x^+r^) 

((x1)2  +  (x2)2  +  (x3+rg)2)3/2 


+  2Q(x®sin  L  -  x®cas  L)  +  Q^x^  (46) 

-  20x‘^sin  L 

+  Q^(x^sin^  L  -  (x^+rQ)cos  L  sin  L)  (47) 
+  2nx^cos  L 

+  (x^+r_)cos^  L  -  x^cos  L  sin  L)  (48) 


where 


p  -  G  m- 


and  all  other  variables  are  defined  in  the  List  of  Symbols.  A  spherical 
harmonic  expansion  of  the  earth’s  gravitational  field  is  probably  not 
necessary  since  most  shots  over  the  range  are  low  altitude  however,  a  major 
effect  will  be  atmospheric  drag  and  more  precise  modeling  may  be  necessary. 

Any  one  of  three  sensor  types  moy  reside  at  each  sensor  location.  The 
three  types  are 

■  range  radar  wherein  range  (r^),  azimuth  (P^)  and  elevation  (9^)  data 
are  available, 


optical  tracker  wherein  azimuth  and  elevation  data  are  available,  and 
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•  doppler  radar  wherein  range,  range  rate  (ri),  azimuth  and  elevation 
data  are  available. 

The  observational  equations  in  terms  of  the  global  state  are 

r^  -  ((x''-d^-l)2  +  (x2-di'2)2  +  (^3_jji,3)2j1/2  ( 


x^-d^’ ^ 

[0  0  1]  x2-di-2 

x^-di • ^ 


0^  »  sin“ 


((x^-d^*^)^  +  (x^-d^*^)^  +  (x^-d^’ 


»  tan“^ 


x^ -d^ ’ ^ 

[0  1  0]  x2-d^-2 

x3_di.3 


X ^ ^ 

[1  0  0]  x2-di-2 

x3_di.3 


xl-di-l 

[  x'^  X®  X®  ]  x2-d^>2 

x^-d^’ 3 


((x^-d^’'')2  +  (x2-di*2)2  +  (x*-di-3)2)1/2 
The  various  measurement  variables  are  defined  by  Figure  1  that  follows. 


,i.1  pi 


Figure  1:  The  Local  Coordinate  System  and  Measurement  Variables 


The  nonlinear  equations  of  motion  (46)-{48)  were  integrated  forwards  in 
time  starting  from  (30)  and  using  the  Fourth  Order  Runge  Kutta  method  with  a 
1  second  interval.  Figures  2  and  3  show  the  resultant  trajectory  of  the 
target  vehicle.  Comparison  with  the  flat  earth  model  predictions  show  close 
agreement  in  peak  altitude  (110,305.  feet),  terminal  time  (166+  seconds),  and 
the  projected  path  length  (82.7  miles). 
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44,403. 

226,391 . 

110,283. 

285.5 

2619.5 

38.7 

83 

44,688. 

229,010. 

110,305. 

285.6 

2619.1 

6.9 

84 

44,975. 

231,629. 

110,296. 

285.8 

2618.7 

-24.9 

166 

69, 187. 

444,565. 

1176. 

308.6 

2570 . 3 

-2640.9 

167 

69,496. 

447,135. 

-1481 . 

309.0 

2569.6 

-2673.0 

Table  2:  Nominal  Positions  and  Velocities  for  Key  Time  Points  of  the 
Ballistic  Trajectory 

As  a  further  check  on  the  Fortran  code  for  this  part  of  the  simulation,  the 
total  translational  energy  (kinetic  plus  potential)  of  the  target  vehicle  was 
computed  for  each  point  along  the  computed  trojectory.  The  total  energy 
remained  constant  to  within  I <  as  it  should  since  there  are  no  external  forces 
(u(t)  »  0)  and  the  system  is  conservative.  Figures  4  through  8  show  some  of 
the  corresponding  measurements  for  the  various  data-types  with  v^  •  0  for 
(UiiM). 

Each  row  of  A(t)  is  computed  by  partial  differentiating  the  seme  row  of 
f(x(t))  w.r.t.  the  nominal  state.  Thus,  A(t)  has  the  following  structure 


A(t)  - 


0 

0 

0 


0 

0 

0 


X 

ay 

ix6 

ix1 

I 

X 

I 

(54) 


where 


Figur«  4:  Sensor  #1  Azimuth  and  Elevation  Measurements  Corresponding  to  the 
Nominal  Ballistic  Trajectory. 
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Figure  5:  Sensor  #2  Azimuth  and  Elevation  Measurements  Corresponding  to  the 

Nominal  Ballistic  Trajectory. 


Figure  6:  Sensor  #3  Range,  Azimuth  and  Elevation  Measurements  Corresponding 

to  the  Nominal  Ballistic  Trajectory. 
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Figure  7: 


Sensor  #4  Azimuth  and  Elevation  Measurements  Corresponding  to  the 
Nominal  Ballistic  Trajectory. 


Figure  8:  Sensor  #5  Range,  Azimuth  and  Elevation  Measurements  Corresponding 

to  the  Nominal  Ballistic  Trajectory. 
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Figure  9:  Perturbed  State  Positions  Versus  Time  for  the  Ballistic  Target. 
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Each  row  of  C^(t)  is  computed  by  partial  differentiating  the  same  row  of 
hi(x(t)  w.r.t.  the  nominal  state.  Thus,  when  the  ith  sensor  is  a  doppler 
radar,  C^(t)  has  the  fallowing  structure 
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where  u  is  the  argument  of  sin“^  in  equation  (51). 
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The  6  biases  a^, d^’ ^ , d^’ d^’ ^  per  local  system  could  bo  included 
as  states  in  the  filter  and  estimated  in  order  to  correct  for  any  preflight 
miscalibration .  Then,  =  0,  =  0  . . .  would  bo  added  to  the  equations  of 

motion  and  additional  partial  derivative  expressions  would  be  needed. 

However,  only  the  6  positions  and  velocities  of  the  target  vehicle  were 
included. 

In  order  to  create  simulated  data  and  test  the  state  estimation  part  of 
the  code,  a  subroutine  RANDOM  for  generating  sequences  of  Gaussian  random 
vectors  with  prescribed  covariance  was  written.  Two  algorithms  were 
considered.  The  first  proceeds  by  rotating  coordinotes  to  o  system  in  which 
the  covariance  matrix  is  diagonal.  In  this  system  the  multivariate  normal 
density  becomes  equal  to  the  product  of  its  marginal  densities,  and  each 
marginal  density  can  be  sampled  independently  of  the  other  components.  After 
obtaining  o  sample  vector  in  this  rotated  system,  the  coordinotes  are  rotated 
back  to  the  original  system. 

The  second  algorithm  proceeds  by  decomposing  the  multivariate  normal 
density  into  the  product  of  the  marginal  density  of  the  first  variate  times 
the  joint  density  of  the  remaining  variates,  conditional  upon  the  value 
sampled  for  the  first.  This  joint  density  is  determined  once  the  first 
variate  has  been  sampled  from  its  marginal  density.  The  procedure  is  then 
applied  to  the  second  variate  and  iterated  until  values  have  been  assigned  to 
all  components  of  the  sample  vector.  This  "Conditional  Decomposition 
Algorithm"  will  execute  more  rapidly  than  the  latter  "Matrix  Diagonalization 
Algorithm"  especially  for  time  varying  covariance  matrices.  Thus  it  was 
chosen  as  the  basis  for  subroutine  RANDOM. 

Assuming  a  constant  covariance  matrix,  RANDOM  was  tested  by  counting  the 
number  of  random  values  within  several  bands  for  each  component.  Comparison 
with  theory  has  shown  agreement  to  within  a  few  percent.  Much  time  was  spent 
in  validating  RANDOM  since  its  accuracy  is  a  prerequisite  for  meaningful 
future  comparisons  of  DSRIF  estimates  with  decoupled  Kalman  filter  estimates 
in  Phase  II  research. 

Figures  9  and  10  show  the  evolution  of  the  perturbed  state  as  governed  by 
equation  (3).  The  initial  condition  is 


Figure  10:  Perturbed  State  Velocities  Versus  Time  for  the  Ballistic  Target 


Figure  11:  RMS  Global  Position  Estimate  Error  Versus  Time. 


Sxq  =  [  o.ft  O.ft  O.ft  2.8629  ft/s  26.4021  ft/s  26.558  ft/s  ] 

which  represents  a  scaling  of  the  initial  velocity  vector  for  the  nominal 
trajectory.  Thus,  the  combined  nominal  plus  perturbed  solution  should  reach 
impact  along  the  nominal  line  of  sight  (which  is  the  vector  difference  between 
the  nominal  impact  point  and  the  nominal  initial  condition)  but  further  from 
Xq  since  a  positive  scaling  factor  was  used.  A  diagonal  with  small 
variances  was  used  to  generate  the  process  noise  sequence.  The  position 
variances  were  .1  ft^  and  the  velocity  variances  were  .1  (ft/s)^. 

Figures  11,  12  and  13  are  the  corresponding  DSRIF  results  with  all  of  the 
prior  and  process  noise  information  embedded  in  the  merge  processor.  Each  of 
the  12  local  filters  processed  1  measurement  variable.  A  diagonal  R^,  with 
variances  of  10”®  deg^  and  10^  ft^  for  angular  variables  and  range  variables 
respectively,  was  used  to  generate  the  measurement  noise  sequence.  The 
initial  state  estimate  for  all  of  the  local  filters  was 

Sxj  »  [20. ft  20. ft  20. ft  23.8629  ft/s  47.4021  ft/s  47.558  ft/s  ] 

and  a  diagonal  Pq(+),  with  variances  of  100  ft^  and  100  (ft/s)^  for  position 
variables  and  velocity  variables  respectively,  was  used  to  initialize  the 
merge  processor.  Figures  11  and  12  show  that  the  rms  position  and  velocity 
estimate  errors  quickly  decay  ta  steady  state  mean  values  after  less  than  10 
time  samples.  The  corresponding  estimate  error  covariances  follow  the  same 
course  as  expected. 

In  Figure  14  the  process  noise  levels  were  multiplied  by  10  and 
comparison  with  Figure  13  shows  that  the  corresponding  estimate  error 
covariances  increase  as  well.  This  is  as  expected  since  Q|^  is  linearly 
related  to  the  time  updated  estimate  error  covariance  i.e.,  the  conventional 
cavariance  time  update  equation  is  given  by 

Pk+l(-)  =  Vk(  +  )*k^''  +  Qk  (93) 

Furthermore,  the  same  phenomenon  results  when  R^  is  multiplied  by  a  factor  of 
5  in  Figure  15.  The  conventional  covariance  measurement  update  equation  in 
Josephson  Stabilized  form 

Pk(  +  )  -  [I  -  P^{-)  [I  -  (94) 

may  be  combined  with  the  Kalman  gain  equation 

Kk  -  Pk(-)Hk^'"  CHkPk(-)H^^'‘  +  R,^]-''  (95) 

to  show  that  the  time  updated  estimate  error  covariance  is  linearly  related  to 
R|^  as  well.  The  result  is  that 

Pk”^(  +  )  -  Pk'^(”)  ^k^'*  '^k“^  ^^k  (96) 


2.3  Extended-Decentralized  Square  Root  Information  Filtering  of  MLRS  Data 


On  November  11,  1987  six  rockets  were  launched  sequentially  in  time  over 
a  period  of  2  hours  and  30  minutes  at  WSMR.  Only  1  rocket  was  airborne  at  any 
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Figure  12:  RMS  Global  Velocity  Estimate  Error  Versus  Time. 
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Figure  13:  Global  Position  and  Velocity  Estimate  Error  Covariances  Versus 
Time.  R^(j.j)-10“  for  angle  variables  and  10^  for  range 
variobles.  Q^(j,j)  -  . 1  f or  all  j-l . 6. 
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Figure  14:  Global  Position  and  Velocity  Estimate  Error  Covariances  Versus 
Time.  Rt<( j  I  j  )“10~®  for  angle  variables  and  10^  for  range 
variables.  Q^(j,j)  -  1.  for  oil  j-1 . 6. 

20  ft  squared 


Figure  15;  Global  Position  and  Velocity  Estimate  Error  Covariances  Versus 

Time.  R|<(j,j)-5  x  10“®  for  angle  variables  and  5  x  10^  for  range 
variables.  O^Cj-j)  -  1.  for  all  j-1 . 6. 


one  time  and  thus  data  association  for  multitarget  tracking  would  not  be 
needed.  MTI  obtained  a  copy  of  the  Multiple  Launch  Rocket  System  (MLRS)  data 
in  order  to  gain  experience  at  processing  real  test  range  single  target  data 
with  the  algorithm.  The  digitized  measurements  for  all  six  shots  were  plotted 
in  order  to  select  the  best  set  as  characterized  by  the  least  amount  of  data 
drop  outs  and  outlyers.  Figures  16  through  26  comprise  the  data  set  chosen 
and  this  was  the  second  shot. 

The  MLRS  data  set  contained  azimuth  and  elevation  angle  measurements 
(with  respect  to  each  local  sensor)  from  11  optical  trackers  located  at  the 
range  coordinates  listed  in  Table  3  below.  The  origin  of  the  WSCS  is 
[500,000.  500,000.  0.]  in  units  of  feet  with  its  latitude  and  longitude 

equal  to  32.38  and  106.481  degrees  respectively.  Each  optical  tracker  is  also 
characterized  by  a  set  of  3  Euler  angles  a^,  and  which  define  its 
orientation  relative  to  the  Global  Coordinate  System  (GCS)  as  in  section  2.2. 
The  first  two  angles  are  calculated  using  equations  (39)  through  (45)  and  the 
angular  misalignment  is  assumed  equal  to  0. 
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-201,753.4857 
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11 
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_ 1 

Table  3:  Optical  Tracker  Locations  and  Orientations  in  GCS  for  the  MLRS  Test 


The  data  set  also  contained  range,  azimuth  and  elevation  angle 
measurements  from  3  radars  but  with  respect  to  the  local  coordinate  system 
originating  at  the  launcher.  Their  locations  and  orientations  are  given  in 
Table  4  below  however,  the  launcher  location  d^'^  =  -14,074.34, 
d^’^  »  -247,569.51,  d^’^  »  2,505.08  was  used  instead. 
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ri 
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-13,971 .87 

-264,430.55 

2,334.06 
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13 
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-27,098.86 
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106.4212 
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14 

394 

-  2,382.95 

-260,481.11 

2,403.07 
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106.3410 

0. 

Table  4:  Radar  Locations  and  Orientations  in  GCS  for  the  MLRS  Test 
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Figure  16;  Rador  #350  Range,  Azimuth  and  Elevation  Measurements  for  MRLS 
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Figure  17:  Radar  #393  Range,  Azimuth  and  Elevation  Measurements  for  MRLS 
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Measurements  for  MRLS. 
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An  extended  version  of  the  DSRIF  is  needed  since  the  nominal  rocket 
trajectories  were  unavailable.  The  Extended  DSRIF  (E-DSRIF)  may  be  derived  by 
extending  the  observations  equation,  linearized  about  the  current  estimate,  to 
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zj  -  b^(x)  -  Xk(-) 
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and 

^  hi(x) 
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and  the  dynomics  equation,  linearized  about  the  current  estimote.  to 

’‘k+1  =  ^k  ’'k  +  ®k  '^k  +  9k 

where 

gk  =  f(x)  -  Fk  Xk(+) 

X«Xk(+) 


(100) 


(101) 
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^  X 


(102) 


|x=Xk(+) 


Substituting  equations  (97)  and  (100)  into  the  least  squares  performance 
functional  of  equation  (14)  gives  the  following  set  of  equations  for  updating 
the  local  filters  in  extended  form: 
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Time  Update 
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Herein,  and  are  orthogonal  transformations  which  put  the  matrices 

on  the  left  hand  sides  of  (103)  and  (104)  into  upper  triangular  form  (they  may 
be  implicitly  computed  using  Householder  transformations).  All  other 
variables  are  defined  in  the  List  of  Symbols. 

Processing  on  the  global  scale  is  the  some  os  for  the  DSRIF  i.e.,  the 
merge  steps  are  exactly  as  defined  in  section  2.1.  Only  processing  on  the 
local  scale  is  modified  as  such.  A  major  difference  between  the  E-DSRIF  and 
DSRIF  is  that  the  local  E-DSRIFilters  require  knowledge  of  the  globally 
optimal  estimate  X(<(+)  in  order  to  compute  their  first  order  Toylor  series 
expansion  terms  Fj^,  g^^,  hJ  and  zj^  whereas  the  DSRIF  may  compute  X|^(  +  )  at  any 
rate  less  than  the  highest  data  rate.  Future  research  should  examine  whether 
an  E-DSRIF  algorithm,  in  which  the  Taylor  series  expansions  are  about  the 
locolly  optimal  estimates,  may  be  derived. 

In  order  to  derive  a  suitable  dynamical  model  as  well  as  initialize  the 
filter,  the  position,  velocity,  acceleration  and  jerk  of  the  rocket  were 
precomputed  using  finite  differencing  with  At  »  . 1  seconds.  Results  are 
plotted  in  Figures  27  through  30  using  all  of  the  dota  provided  for  radar  #350 
except  for  the  first  21  samples  (we  estimated  that  all  rt’s  hod  locked  onto 
the  target  by  the  22nd  sample).  Figure  30  indicates  that  jerk  is  suitably 
modeled  as  a  white  Gaussia  noise  process  with  constant  mean.  Thus,  the  E- 
DSRIF  was  encoded  in  Fortran  '77  using  a  second  order  polynomial  dynamical 
model  wherein 
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from  Radar  #350. 


eration  Measurements  from  Radar  #350. 
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The  sample  mean  and  covariace  of  was  computed  to  be 

[  -10.  ft/sec^  -.055  ft/sec^  -.48  ft/sec^  ] 

and 

diag  [  23,851.  ft^/sec®  1,458.  ft^/sec®  23,199.  ft^/sec®  ] 

respectively,  however,  a  much  larger  covariance  was  used  in  order  to  compensate 
for  any  errors  in  the  model.  Detailed  testing  of  the  algorithm  is  deferred  to 
Phase  II  work  wherein  an  adaptive  method  for  adjusting  0^  in  real  time  will  be 
investigated . 

Figures  31  through  35  show  the  rt  and  ot  measurements  as  predicted  by  the 
E-DSRIF.  Comparison  with  the  actual  measurements  in  Figures  16,  18,  19,  20 
and  22  shows  on  exact  match  to  within  a  plotting  line  width.  A  better  means 
of  comparison  is  thus  provided  below  in  Tables  5  and  6. 


Figure  32:  Radar  #394  Range,  Azimuth  and  Elevation  Estimates  for  MLRS. 

R|<(  j  I  j  for  ot  variables  and  for  rt  range,  azimuth 

and  elevation  variables  respectively. 

Qk(j.j)  -  diag  [  1.7  x  lO"''  7.4  x  lO''^  9.3  x  lo''^  ]  . 

590/58 


X. 


i — 


X. 


58sec 


Figure  33:  630  Azimuth  and  Elevation  Estimates  for  MLRS. 

Rk(  j .  j  ^or  ot  variables  and  10.,1.,1.  for  rt  range,  azimuth 

and  elevation  variables  respectively. 
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11.85 

459  m  80,328.90  359.22 

e  80,254.10  359.22 


80, 178.71 
80,254.10 


359.22 

359.22 


Table  5:  Radar  #350  and  #394  Measurements  and  Estimated  Measurements  for 

MLRS.  R(<(  j .  j  )*10^°  for  ot  variables  and  10.,1.,1.  for  rt  range, 
azimuth  and  elevation  variables  respectively. 

OkCj.j)  =  diag  [  1.7  x  lO'''^  7.4  x  lo"'^  9.3  x  lO'''^  ]  . 
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255.49 
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m 

154.39 

1.19 

150.73 

3.58 

258.16 

4.07 

e 

154.90 

.25 

151.71 

2.14 

256.30 

3.00 

229 

m 

68.57 

28.54 

22.78 

21.75 

««« 

««« 

e 

69.20 

26.13 

22.93 

20.28 

331.64 

12.52 

230 

m 

«  «  « 

««« 

22.64 

21 .65 

««« 

««« 

e 

68.74 

26.07 

22.80 

20.18 

331 .73 

12.48 

459 

m 

26.07 

5.58 

#«« 

«»« 

«  «  « 

e 

26.32 

3.94 

10.98 

3.66 

342.28 

2.98 

_ 

Table  6:  G30,  G80,  G110  Measurements  and  Estimated  Measurements  for  MLRS. 

variables  ond  for  rt  range,  azimuth 

and  elevation  variables  respectively. 

Qk(J.j)  =  'liag  [  X  lO"'^  7.4  X  lO''^  9.3  x  lo''^  ]  . 

***  denotes  data  drop-out 

The  large  values  of  Rk(j,j)  for  ot  variables  serves  to  weight  the  rt  data  much 
more  heavily  in  computing  estimates.  Decreasing  the  ot  measurement  errors  to 
more  realistic  values  should  give  similar  results  since  the  predicted  ot 
measurements  matches  their  actual  values  very  closely. 

Figures  36  and  37  show  the  global  position  estimates  and  corresponding 
estimate  error  covariances  respectively.  Agoin,  the  rocket  positions  derived 
from  radar  #350  as  compared  with  the  E-DSRIF  estimates  based  upon  all  of  the  5 
selected  sensors,  show  extremely  close  agreement.  The  slight  difference  in 
the  estimate  of  height  is  due  to  using  »  I  instead  of  its  correct 

value  as  defined  by  equation  (34).  In  Table  7  below,  the  estimates  are 
compared  using  the  correct  cooordinate  tronsformation . 
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459  e  -15,197.  -167,388.  5,729. 

m  -15,202.  -167,307.  5,704. 


-431.  21.  -221.  184. 

-439.  .9  -25.  59. 


Table  7: 


Global  Position  Estimates  and  Derived  Measurements  for  MLRS. 

j  •  j  )  =  ot  variables  and  10.,1.,1.  for  rt  range,  azimuth 

and  elevation  variables  respectively. 

Ou(j.j)  =  diag  [  1.7  x  lO’’'  7.4  x  lo’’’^  9.3  x  lo'’’^  ]  . 


Finally,  the  monotonically  increasing  estimate  error  covariance  (  actually. 


Pq  =  diag  [  10^^  ft^  ...  10^^  ft^/sec^  ...  10^^  ft^/sec*^  ...  ]  was  used  to 

initialize  the  covariance  propagation  so  that  the  first  step  is  a  large,  off 
scale  decrease  to  approximately  100  ft^)  is  due  to  our  using  values  of  Qj^ 
approximately  12  orders  of  magnitude  higher  than  its  precomputed  sample  value . 
A  more  realistic  value  should  result  in  a  P|^  with  quite  the  opposite  behavior. 


was  used  to 


3.0  Estimates  of  Technical  Feasibility 

The  objective  of  Phase  I  research  was  to  determine  the  feasibility  of 
constructing  an  integrated  test  range  tracking  network  based  upon  the  DSRIF. 

A  multitude  of  test  range  scenarios  is  envisioned  so  thot  a  robust  system  is 
needed.  At  one  extreme,  test  vehicles  may  include  ballistic  projectiles  with 
well  defined  nominal  trajectories  a  priori  while  ot  the  other,  multiple  smart 
munitions  with  maneuvering  capability  is  possible.  The  key  to  a  successful 
network  design  is  to  employ  a  more  or  less  sophisticated  version  of  the 
algorithm  depending  upon  the  particular  scenario.  Thus  the  network  must  be 
adaptable.  For  example,  preflight  simulations  of  the  proposed  shot  using  high 


fidelity  aerodynamic  models  can  yield  good  values  for  the  process  noise  levels 
and  a  basic  DSRIF  should  result  in  good  tracking  performance.  However,  a 
sudden  departure  from  the  nominal  trajectory  would  require  a  detection 
mechanism  as  part  of  the  algorithm  and  adjustment  of  in  real  time. 

In  order  to  determine  the  feasibility  of  our  distributed  approach  to 
multisensor  tracking,  several  specific  technical  objectives  must  be  met. 

First  and  foremost,  the  basic  DSRIF  theory  needs  to  be  extended  to  enable  the 
tracking  of  maneuvering  vehicles,  high  dynamic  trajectories,  and  multiple 
targets.  The  latter  requires  that  a  theory  for  associating  data  with  targets, 
based  upon  the  DSRIF,  be  developed.  Other  theoretical  questions  such  as  the 
development  of  a  delayed-state  DSRIF  for  processing  range-rate  measurements,  a 
method  for  isolating  faulty  sensors,  and  efficient  implementations  of  the 
DSRIF  that  facilitate  high  data  rates  need  to  be  addressed. 

Secondly,  the  DSRIF  is  a  new  algorithm  which  has  undergone  only  limited 
testing  in  Phase  I  research.  Extensive  testing  within  a  multisensor 
multitarget  tracking  scenario  is  needed.  Finally,  consideration  needs  to  be 
given  to  the  design  of  the  tracking  network  both  at  the  global  and  local 
levels.  The  major  question  here  is  whether  a  sufficient  data  rate  can  be 
achieved  using  current  chip  technology.  Another  question  is  whether  the 
architecture  can  be  reconfigured  (in  software)  to  implement  other  members  of 
the  family  of  DSRIFs.  A  DSRIF  based  multisensor  laboratory  tracking 
experiment  should  be  performed. 


For  multitarget  tracking,  correlation  of  meosurements  with  targets  con 
best  be  done  using  a  hypothesis  testing  approach.  The  idea  is  to  select  the 
correlation  of  measurements  with  targets  thot  has  maximum  probability  given 
the  data.  Calculation  of  all  combinations  to  form  the  entire  set  of  these 
conditional  probabilities  can  be  prohibiting,  especiolly  in  a  dense  target 
environment.  A  major  advantage  in  using  the  DSRIF  is  the  tremendous  reduction 
in  computational  cost  associated  with  this  calculation. 
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